################################################################################
# Stat-Tutorial-05-LawLargeNumbers.r                                           #
#                                                                              #
# This is a tutorial on the law of large numbers.                              #
#                                                                              #
# R. Labouriau                                                                 #
# Last revised: Fall 2020                                                     #
################################################################################
# Copyright © 2019 by Rodrigo Labouriau

################################################################################
# The law of large numbers for a binomial distributed variables
################################################################################

# Here we set the simulation parameters.
n.rep <- 1000              # Number of simulations
n.trials <- 20             # Number of binomial trials (or size)
prob.success <- 1/2        # Parameter describing the probability of success in 
                           # each binary trial
x <- rbinom(n=n.rep, size=n.trials, prob=prob.success)
x
mean(x)
n.trials*prob.success
mean(x) - (n.trials*prob.success)

# Studying the behaviour of the mean for an increasing number of replications
X <- numeric(1000)
N.rep <- 1:1000 + 3 ; N.rep
for(i in 1:1000){
  bb <- rbinom(n=N.rep[i], size=n.trials, prob=prob.success)
  X[i] <- mean(bb)
}
# pdf("~/Dropbox/Courses/12-StatisticalModelling/Lectures/FigureLectureDay2-01.pdf")
plot(N.rep, X, ylab='Sample mean', xlab='Number of repetitions')
abline(h= n.trials*prob.success, col='red', lwd=3)
# dev.off()

################################################################################
# Studying the behaviour of the mean for an increasing number of replications
# for the Poisson
################################################################################

theoretical.mean <- 10
X <- numeric(1000)
N.rep <- 1:1000 + 3 ; N.rep
for(i in 1:1000){
  bb <- rpois(n=N.rep[i], lambda=theoretical.mean)
  X[i] <- mean(bb)
}
plot(N.rep, X, ylab='Sample mean', xlab='Number of repetitions', main='Poisson')
abline(h=theoretical.mean, col='red', lwd=2)

for(i in 1:1000){
  bb <- rpois(n=N.rep[i], lambda=theoretical.mean)
  X[i] <- mean(bb)
}
points(N.rep, X)

for(kk in 1:10){
  for(i in 1:1000){
    bb <- rpois(n=N.rep[i], lambda=theoretical.mean)
    X[i] <- mean(bb)
  }
  points(N.rep, X)
}

# Studying the behaviour of the mean for an increasing number of replications
theoretical.mean <- 10
X <- numeric(1000)
N.rep <- 1:1000 + 3 ; N.rep
for(i in 1:1000){
  bb <- rpois(n=N.rep[i], lambda=theoretical.mean)
  X[i] <- mean(bb)
}
plot(N.rep, X, ylab='Sample mean', xlab='Number of repetitions')
abline(h= n.trials*prob.success, col='red', lwd=2)

################################################################################
# Studying the behaviour of the mean for an increasing number of replications
# for the Normal distribution.
################################################################################

theoretical.mean <- 0
X <- numeric(1000)
N.rep <- 1:1000 + 3 ; N.rep
for(i in 1:1000){
  bb <- rnorm(n=N.rep[i], mean=theoretical.mean)
  X[i] <- mean(bb)
}
plot(N.rep, X, ylab='Sample mean', xlab='Number of repetitions', main='Normal')
abline(h=theoretical.mean, col='red', lwd=2)

for(kk in 1:10){
  for(i in 1:1000){
    bb <-  rnorm(n=N.rep[i], mean=theoretical.mean)
    X[i] <- mean(bb)
  }
  points(N.rep, X)
}
abline(h=theoretical.mean, col='red', lwd=3)

# THE END !